当红辣子鸡-空间转录组与单细胞转录组的整合分析（上篇）

空间转录组技术发展史

如果说最近什么最红火，那一定是奥运会，今年的奥运会从开始到现在就一直吸引着国人的注意，奥运会的精彩让我差点忘记了地球外太空还有我们的国宝-“三位航天员”；在科研界，当红辣子鸡是谁呢？那肯定属于空间转录组，空间转录组被Nature Methods杂志评为2020年度技术，今年空间转录组技术得到的广泛的应用，已经成为了当前科研界名副其实的“当红辣子鸡“。

何为空间转录组？空间转录组与单细胞转录组有何区别？空间转录组学在空间上进行转录组分析，这保留了空间位置的转录组学，让我们对实体组织空间上的异质性进入了一个全新的视角，单细胞组学是细胞分辨率上让我们了解组织样品的细胞异质性，可以让我们知道样品中的细胞图谱。被Science评为2018年年度突破技术的单细胞组学技术和被Nature Methods评为2019年年度技术方法的单细胞多组学技术让我们在了解生物异质性和从更高分辨率较多解决生物学问题提供了极大的方便，但是令科学界遗憾的是，这种单细胞组学技术没有空间位置信息，就像我们知道银河系有哪些行星，但是如果我们不知道其具体的位置信息的话，这对科学界了解事物本质问题还是有一定的困扰。因此空间转录组技术的诞生给我们提供了从空间位置了解事物本质方法途径。目前空间转录组技术基本上还未到单细胞水平，一般一个空间位置spots上有多个细胞，10X Genomics Visum一般为1-10个细胞，空间转录组学能够提供空间位置，但是分辨率达不到单细胞水平，单细胞转录组学分辨率能够达到单细胞分辨率，但是没有空间位置信息，因此将空间转录组与单细胞转录组结合的话，那就既有空间位置信息，又达到了单细胞分辨率，这为科学界提供了有力的工具，大大的提高了生物体的组织、器官等研究的分辨率和准确性。

其实空间转录组学的发展历史可以追溯到上个世纪，经历过漫长的发展，从早期的smFISH（单分子RNA荧光原位杂交技术）和LCM（激光线为切割）技术发展到目前的10X Visium（空间分辨的全转录组信息）、Slide-Seq（溶液中的珠子）、APEX-Seq（亚细胞定位），从早期的少量基因到目前的全转录组（所有的基因），空间转录组技术在近30年中得到的极大的发展，特别是随着10X Genomics公司收购Visium以及后续10X Visium产品的推出，这极大的推动了空间转录组技术的应用。

10X Visium空间转录组基本原理

10X Visium空间转录组是目前应用最广的空间转录组技术，其基本原理如下：一般有两种载玻片，一种载玻片是组织透化条件优化玻片，另外一种载玻片为文库构建载玻片。组织透化条件优化载玻片是主要在建库之前，不同组织和样品的透化条件可能不一样，需要根据自己的样品组织进行透化条件的优化，以达到最佳透化条件进行后续建库，这里主要透化优化包括透化时间和透化浓度等优化。在建库载玻片上有四个捕获区域，每个捕获区域的大小为6.5x6.5 mm，包含5000个spots（被条形码标记的点），每个spots的直径为55 μm，spots和spots之间中心的距离为100 μm，并且每个spots都有一个唯一的barcode序列，这与单细胞的细胞barcode一样。一般将切割好的样品平铺在载玻片中间位置，然后根据优化获得透化条件进行透化，组织样品中细胞会释放mRNA，然后与每个spots的polyT序列进行结合，然后进行后续的文库构建测序，根据spots的barcode确定空位位置，从而实现空间转录组学的检测。需要注意的是，这里的spots并不是单细胞级别的，一般含有1-10个细胞。

10X Visium空间转录组基本原理

空间转录组与单细胞转录组整合分析

器官系统由不同的细胞亚群组成，它们在给定组织内的空间位置与其功能密切相关。单细胞转录组表征单个细胞的转录组，可以揭示给定器官内的细胞亚群。然而，在单细胞转录组必要的组织解离步骤中，单个细胞的分离破坏了它们在天然组织中的空间定位及其彼此接近的信息。鉴于近分泌和旁分泌信号的作用范围为 0 到 200 μm，因此此类空间信息对于了解正常组织和患病组织的细胞间通讯至关重要。使用空间转录组学检测完整组织，通过物理定位在单细胞转录组识别的特定细胞亚群中表达的基因集解决了这一挑战。当前的空间转录组学方法本身还不能提供组织中精确定位的单细胞的深层转录组信息。然而，它们可以阐明富含不同基因集的空间位置。这将细胞亚群定位到它们所在的组织“邻域”，并指定它们表达的配体和受体以影响局部细胞间通讯。因此，当结合使用时，单细胞转录组和空间转录组学可以在其天然组织环境中定位具有转录特征的单个细胞。因此，整合单细胞转录组和空间转录组学数据可能会增加我们对特定细胞亚群的作用及其在发育、稳态和疾病中的相互作用的理解。

空间转录组与单细胞转录组整合分析的方法随着空间转录组技术的发展，已经有了很多种方法，目前主要有两大类：

1. 映射方法：映射方法就是将单细胞转录组 细胞映射到用户定义的区域或2D或3D空间中的单细胞坐标。然而，这些方法依赖于具有原型结构的组织。

2. 解卷积：就是认为空间转录组是spots是由混合细胞组织，通过解卷积方法将spots是细胞组成给预测出来，然后根据细胞类型的比例来进行后续研究。

其实空间转录组与单细胞转录组整合分析，在空间转录组分析上就是空间转录组上细胞亚群鉴定的一个过程，细胞亚群的鉴定是单细胞转录组分析中最为关键的一个步骤，同样的空间转录组上，其细胞亚群的鉴定（空间转录组与单细胞转录组整合分析）也是空间转录组分析的关键步骤，也是下游的细胞通讯、发育轨迹分析的基础。

空间转录组与单细胞转录组整合方法的映射方法，其实很多属于多组学的整合方法，已经有很多成熟的工具用于其整合分析，比如LIGER、Seurat Integration（ Seurat 3.0以上）和Harmony，这些整合方法得到的结果大部分是一个spots属于某个细胞类型，但是实际上，每个spots是属于多细胞，一般含有不同类型的细胞，因此这种映射方法在很多分析中应用不是很好。基于混合细胞的假设的解卷积的方法就很好的解决了这个问题，这种方法更符合研究实际，因此在这里，主要介绍解卷积的方法，这里介绍两种工具：Robust cell-type decomposition （RCTD）和SPOTlight。

RCTD

RCTD（Robust cell-type decomposition），是一种将 RNA 测序混合物分解为单个细胞类型的监督学习方法，从而能够将细胞类型分配给空间转录组像素。具体来说，我们利用带注释的 scRNA-seq 数据来定义空间转录组学数据中预期存在的细胞类型的细胞类型特异性概况。

有监督的细胞类型学习的一个相关挑战是常说的平台效应，即依赖于技术的文库制备对测序平台之间单个基因捕获率的影响。研究表明，如果不考虑这些平台效应，监督方法不太可能成功，因为系统技术可变性主导了相关的生物信号。这些影响先前已在同一生物样本上的单细胞和单核 RNA-seq (snRNA-seq) 之间的比较中发现，但是RCTD解决了其平台效应，可以应用于不同平台的空间转录组的解卷积分析。

RCTD 可以准确地发现模拟和真实空间转录组数据中细胞类型的定位。此外，我们表明 RCTD 可以检测细微的转录组差异以在空间上绘制细胞亚型。最后，我们使用 RCTD 来计算预期的细胞类型特异性基因表达，从而能够根据细胞的空间环境检测基因表达的变化。

软件的安装

最近安装GitHub上的包老是不稳定，如果这样安装报错的话，最好还是下载RCTD（https://github.com/dmcable/RCTD）到本地以后再安装吧。

1# install.packages("devtools")
2devtools::install_github("dmcable/RCTD", build_vignettes = TRUE)

软件测试

包的导入

1library(RCTD)
2library(Matrix)
3library(Seurat)
4library(ggplot2)

数据准备

单细胞转录组数据的准备，一般是用Seurat4.0（或者Seurat3.0以上版本）分析的单细胞转录组数据的对象，保存为RDS格式。通常单细胞对象中含有已经鉴定好的细胞类型的结果，保存在[e-mail]矩阵中，列名为cell_annotation，当然也可以改成其他名称。最后生成一个参考数据对象，并保存为RDS格式，便于下次使用和中间文件的保存。

 1sc<-readRDS('sc.celltype.rds') 
 2### Load in/preprocess your data, this might vary based on your file type
 3#refdir <- system.file("extdata",'Reference/Vignette',package = 'RCTD') # directory for the reference
 4counts <- as.matrix(sc@assays$RNA@counts )#read.csv(file.path(refdir,"dge.csv")) # load in counts matrix
 5#rownames(counts) <- counts[,1]; counts[,1] <- NULL # Move first column to rownames
 6meta_data <- [e-mail] #read.csv(file.path(refdir,"meta_data.csv")) # load in meta_data (barcodes, clusters, and nUMI)
 7cell_types <- meta_data$cell_annotation; names(cell_types) <- rownames(meta_data)#meta_data$barcode # create cell_types named list
 8cell_types <- as.factor(cell_types) # convert to factor data type
 9nUMI <- meta_data$nCount_RNA; names(nUMI) <- rownames(meta_data)#meta_data$barcode # create nUMI named list
10### Create the Reference object
11reference <- Reference(counts, cell_types, nUMI)
12#> Warning in Reference(counts, cell_types, nUMI): Reference: nUMI does not match colSums of counts. If this is unintended, please correct this discrepancy. If this
13#>             is intended, there is no problem.
14## Examine reference object (optional)
15print(dim(reference@counts)) #observe Digital Gene Expression matrix
16#> [1] 384 475
17table(reference@cell_types) 
18saveRDS(reference,'ref.rds')

空间转录组数据的准备，空间转录组数据，一般使用Seurat包进行分析的空间转录组结果对象，如果不是，则需要主要对象中含有空间坐标信息，这是必须的，如果没有，可以通过一定的方式手动添加。这里主要用的数据与单细胞数据一样，是原始表达矩阵数据counts。

 1rds_file='spatial.rds'
 2spatial <- readRDS(rds_file)
 3#https://github.com/dmcable/RCTD/issues/26
 4coords <- spatial@tools$[e-mail][,1:2]#GetTissueCoordinates(spatial, scale = NULL) # pulling unscaled tissue coordinates
 5#https://support.10xgenomics.com/spatial-gene-expression/software/pipelines/latest/output/images
 6#spatialcounts <- as.matrix(GetAssayData(spatial, assay = "Spatial", slot = "counts"))
 7spatialcounts <- as.matrix(spatial@assays$RNA@counts)
 8nUMI <- [e-mail]$nCount_RNA
 9names(nUMI)<-rownames([e-mail])
10### Create SpatialRNA object
11puck <- SpatialRNA(coords, spatialcounts, nUMI)

结果预测

1myRCTD <- create.RCTD(puck, reference, max_cores = 8)
2myRCTD <- run.RCTD(myRCTD, doublet_mode = 'doublet') #doublet_mode = 'doublet'

经过前面数据的准备，这里的结果预测还是比较方便，一般只有两行命令，首先创建RCTD对象，使用函数为：create.RCTD，一般的参数不用调整，这里主要参数为max_cores，这里为线程数目，可以根据自己的空间转录组数据和集群资源进行调整。具体其他参数如下，可以根据自己的需求进行调整。

create.RCTD(
   spatialRNA,
   reference,
   max_cores = 8,
   test_mode = FALSE,
   gene_cutoff = 0.000125,
   fc_cutoff = 0.5,
   gene_cutoff_reg = 2e-04,
   fc_cutoff_reg = 0.75,
   UMI_min = 100,
   UMI_max = 2e+07,
   UMI_min_sigma = 300,
   class_df = NULL,
   CELL_MIN_INSTANCE = 25,
   cell_type_names = NULL
 )

然后再直接运行run.RCTD，这里主要参数为doublet_mode，为模型参数，模型参数一般为doublet，双细胞模型，但是并不是所有的spots都为双细胞模型，因此可以选择其他的模型，比如multi和full，不过一般选择doublet为多，multi表示spots可以有多个细胞类型，full表示可以有任何数目的细胞类型，如果选择multi，则是通过一种贪婪算法，知道最后数目固定。具体说明如下：

If in doublet mode, fits at most two cell types per pixel. It classifies each pixel as 'singlet' or 'doublet' and searches for the cell types on the pixel. If in full mode, can fit any number of cell types on each pixel. In multi mode, cell types are added using a greedy algorithm, up to a fixed number.

结果展示说明

 1results <- myRCTD@results
 2# normalize the cell type proportions to sum to 1.
 3norm_weights = sweep(results$weights, 1, rowSums(results$weights), '/') 
 4cell_type_names <- myRCTD@cell_type_info$info[[2]] #list of cell type names
 5spatialRNA <- myRCTD@spatialRNA
 6resultsdir <- 'RCTD_Plots' ## you may change this to a more accessible directory on your computer.
 7dir.create(resultsdir)
 8# make the plots 
 9# Plots the confident weights for each cell type as in full_mode (saved as 
10# 'results/cell_type_weights_unthreshold.pdf')
11plot_weights(cell_type_names, spatialRNA, resultsdir, norm_weights) 
12# Plots all weights for each cell type as in full_mode. (saved as 
13# 'results/cell_type_weights.pdf')
14plot_weights_unthreshold(cell_type_names, spatialRNA, resultsdir, norm_weights) 
15# Plots the weights for each cell type as in doublet_mode. (saved as 
16# 'results/cell_type_weights_doublets.pdf')
17plot_weights_doublet(cell_type_names, spatialRNA, resultsdir, results$weights_doublet, 
18                     results$results_df) 

这里将结果提取出来，然后进行展示，这里抓手结果主要有三种，分别为：cell_type_weights_unthreshold.pdf、cell_type_weights.pdf、cell_type_weights_doublets.pdf，结果会直接在输出目录下，这里的输出目录为：RCTD_Plots。

cell_type_weights.pdf文件主要是单个细胞类型，分别画图，每个spots代表含有当前细胞类型的所有spots，颜色的深浅代表细胞的占比，越红表示越高，越蓝代表越低。

cell_type_weights_doublets.pdf表示双细胞模型下，某个细胞类型所有spots的比例，颜色的深浅代表细胞的占比，越红表示越高，越蓝代表越低。前面两种结果，都是有一定的阈值，过了阈值以后的结果，当然还有没有阈值的原始结果cell_type_weights_unthreshold.pdf，其颜色含义与前面一致，这里是不管值多低，结果都进行展示。

上面展示了在空间位置上的spots细胞类型结果，接下来展示了空间聚类结果上面每种细胞类型的数目情况，如下图所示，X轴代表每个空间聚类，纵坐标y轴代表属于该颜色细胞类型的spots数目。具体代码如下：

1# Plots the number of confident pixels of each cell type in 'full_mode'. (saved as 
2# 'results/cell_type_occur.pdf')
3plot_cond_occur(cell_type_names, resultsdir, norm_weights, spatialRNA)

这里RCTD结果展示完了，其实RCTD文件结果是一个矩阵，矩阵包含了每个细胞属于所有细胞类型的一个概率值，其和加起来为1，具体结果如下图所示：

需要注意的是，如果单细胞数据细胞类型含有未知的细胞类型，一定要删掉，不然会影响结果。

未完待续

请大家继续关注我们的下篇，对SPOTlight的介绍。

参考文献：
[1] Longo S K ,  Guo M G ,  Ji A L , et al. Integrating single-cell and spatial transcriptomics to elucidate intercellular tissue dynamics[J]. Nature Reviews Genetics.
[2] V Marx. Method of the Year: spatially resolved transcriptomics[J]. Nature Methods, 2021, 18(1):9-14.
[3] Asp M , Bergenstrhle J , Lundeberg J . Spatially Resolved Transcriptomes—Next Generation Tools for Tissue Exploration[J]. BioEssays, 2020, 42(10).
[4] Cable D M , Murray E , Zou L S , et al. Robust decomposition of cell type mixtures in spatial transcriptomics[J]. Nature Biotechnology, 2021:1-10.
[5] Dong R ,  Yuan G C . SpatialDWLS: accurate deconvolution of spatial transcriptomic data[J]. Genome Biology, 2021.

作者

尧小飞

审稿

童蒙

编辑

amethyst

往期回顾